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We study the emergence of pattern formation and chaotic dynamics in the one-dimensional 
(ID) generalized Kuramoto-Sivashinsky (gKS) equation by means of a time-series analysis, in 
particular a nonlinear forecasting method which is based on concepts from chaos theory and 
appropriate statistical methods. We analyze two types of temporal signals, a local one and a 
global one, finding in both cases that the dynamical state of the gKS solution undergoes a 
transition from high-dimensional chaos to periodic pulsed oscillations through low-dimensional 
deterministic chaos with increasing the control parameter of the system. Our results demonstrate 
that the proposed nonlinear forecasting methodology allows to elucidate the dynamics of the 
system in terms of its predictability properties. 
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1. Introduction 


The emergence of complex spatiotemporal behavior in spatially extended systems (SES) as a result of 
several mechanisms n onlinearly interacting with each other has a t tracted a lot of attent ion over the last 


few decades (see e.g. Cross fc Hohenbere . Il993l : iKnoblochl. I 2 OO 8 I : iHouehton et al\ . l201Cll |i. A well-known 


example of SES is the Kuramoto-Sivashinsky (KS) equation: 

du du d'^u d^u . , 

where v \s a, viscosity damping parameter. This equation is one of the simplest prototypes exhibiting 
spatiotemporal chaos that retains some of the basic ingredients of any nonlinear SES, namely instabil¬ 
ity/energy production {d'^u/dx^), stability/energy dissipation [d^u/dx^), and nonlinearity (udu/dx) that 
transfers energy from low to high wave numbers. 

E quation JEL was firs t deri ved in the contexts of phase t urbulence in reac tion-diffusion sys¬ 


tems Kuramoto fc Tsuzukil . 197(^ . wrinkled flame front propagation SivashinskvI . 1977| | and unstable drift 
waves driven by electron collision in a tokamak LaQuev et aL . 1975| |. Subsequently, it has been shown 
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to be applicable in a wide spect r um o f physical settin g s, inc l uding hydrodynamic ( e.g. th in-film) instabil¬ 
ities [Shivashinskv fc Michelsoi] . 198d : Babchin et al\. 1983 : Hooper fc Grimshaw . Il985 | and optics such 
as bright spots formed by self-forcing of the beam profile Munkel fc KaiseiT 19961. For thin film flows in 
particular, such as falling films, the KS equation is obtained yia a weakly nonlinear expansion of the ID 
Nayier-Stokes equations s ubject to wall and free-surface boundary conditions and assuming strong surface 


tension effects-long waves [Homsv . 1974| | (derived before Kuramoto’s paper). Not surprisingly, understand¬ 


ing the dynamical behavior of the KS solution has been a major topic of research. For example, it has been 
shown that decreasing the viscosity damping coefficient u in small system sizes with periodic boundary 
conditions leads to a transition from ste ady states to chaos through a p eriod-doubling bifurcation process, 
which is known as Feigenbaum scenario Smvrlis &: Papageorgiou . 199 ll ]. 

With the addition of dispersion, the KS equation becomes the gKS equation: 


du du d‘^u d^u d'^u 


(2) 


where (5 is a positive parameter that characterizes the relative importance of dispersion. This equa¬ 
tion then embraces all elements of any nonlinear process that involves wave e volution in ID, and has 
been reported in a wide variety of contexts, includi ng a reactive fall i ng fil m 


Trevelyan &: Kalliadasid. 


2nn4l |. a film falling down a uniformly heated wall (Kalliadasis et oI\ . l2r)03l | an d a film falling down a 


vertical fiber Ruver-Quil &: Kalliadasia . 2012| | (see also the recent monograph Kalliadasis et al . 2012 1 


for t h e derivation of the KS a nd gKS equations for th in-film flows), in step dynamics Sato fc Uwaha 


1994 : Misbah &: Pierre-Louis . 19961 : Sato et 1998l |. An extended 
local term has also been der i ved f or conducting liquid film in the presence 
field 


Tselniko fc Papageorgioi]! . l2nid | 


version that includes a non- 
of an external electric 


It is now well-known that for small values of 6 the dynamics of the gKS equation resembles the spa- 
tiotemporal chaotic behavior of the KS solution, while sufficien tly large va l ues o f the dispersion parameter 
tend to arrest it in favor of spatially periodic travelling waves (KawaLral . Fl^ . In a regime of moderate 
values of 5 however, traveling waves or pulses appear to be randomly interacting with each other giving 


rise to what is widely k nown as weak/dissipatiye turbulence (in the “Manneyille sense” Mannevillel . 1985 : 


e.g. coherent-structure theories Balmforthl.ll995l:IChang et aHll995: 

Ei &: Ohta. 1994: Duprat et al.. 2009: 

Tseluiko et al. 2010bllal: iDuorat et al. 

. 2011: Tseluiko & Kalliadasis 

. 20141 which have been extended to 

more involved equations Pradas et al. 

. 2011a. 2012a 1 while the regularizing effects of dispersion have been 

considered bv Chang et al. Chang et al. 

1993 ( 1 . who constructed bifurcation diagrams for the periodic and 


solitary solutions of the gKS equation. 

Hence, the gKS equation constitutes a paradigmatic model for the transition from chaotic to regular 
behaviour in spatially extended systems and as such it has been extensively analysed from a spatiotemporal 
analysis point of view (such transitions are nowadays relatively well-known) but it has not received attention 
from the point of view of pure temporal signal analysis. In fact, there is not up to date a detailed analysis of 
the chaotic properties of the gKS solution based on time-series tools (chaos theory), in particular nonlinear 
forecasting methodologies, and how dispersion affects e.g. its predictability propertie s (there are, on the 


other hand, seve r al stu dies on the chaotic behavior of bot h the KS equation [e.g. Chian et all 12002 
Rempel &: Chian . 2003( | and a damped version of it [see e.g. Elder et al. 


a quation le.g. lunian et al.\. I 2 UU 2 I : 
, 1 19971 : Rempel et al. , 2007| | ). We 


note that temporal signals are very often the only experimentally accessible data, hence the importance of 
understanding their dynamics and behaviour. 

Chaos theory provides a comprehensive and systematic treatment of random-like complex dynamics 
in nonlinear systems arising from a broad spectrum of disciplines. Generally, a chaotic signal is defined in 
terms of the properties of exhibiting short-term predictability and long-term unpredictability, something 
which is a consequence of the strong dependence on the initial conditions. Therefore, given a temporal signal 
which describes a quantity of interest in the system, methodologies like nonlinear forecasting can be used 
to quantify many important invariants, such as the fractal dimension, Lyapunov exponents, and entropies 
for irregular temporal evolutions yielding a physical d escription of the corresponding dynamical structures 
and predictability properties of the system (see e.g. [Casdagli . 19891 : Abarbanel et al . 1993 : Baek et al 
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2004 ]), with many importan t applications in e.g . weather prediction jPatil et al\ . 1200 ll : ISlingo fc Palmed . 


2011 ] and transistor circuits Hanias et al. . l2010l | 


Likewise, another related qu estion of interest i s whether the particular time series can be distinguishe d 
from being stochastic or chaotic [Sugihara fc Mav . 1990l : Cazelles fc Ferriere . 199 Tsonis fc Elsnei . 1992ll . 
We note that stochastic effects are ubiquitous in many systems fsee e.g. [Lauritsen et al . , 1996 ; Nicoli et al. 


20101 : Pradas et al. . 201 Ibl . 2012bl | for a noisy version of the KS equation) and they may be present in any 


low-dimensional description of a nonlinear system, as the recent study on stochastic mode reduction of 
the KS e quation using concepts from evolutionary renormalization group theory and information entropy 
indicates Schmuck et al . 2013al li3|. It is actually in this context that nonlinear forecasting methodologies 
become a valuable tool as the more traditional time-series analysis techniques, such as computing the 
power spectral density (PSD) or looking at the trajectories in the phase space, cannot distinguish between 
the two; it is noteworthy that both high-dimensional chaos and stochastic processes produce completely 
scattered plots in the phase space. Hence, providing a reliable and systematic techniqu e that can answer 


these questions has received considera ble attention over the years Guckenheimer . 1982 : Gao et al . 20061 : 


Rosso et all 120071 : IZunino et ali . l2012l | . 


In this work, we undertake a detailed study of the characteristics of the dynamical state of the gKS 
equation in a wide range of the parameter 6 by making use of a time-series analysis based on chaos theory, 
and in particular a nonlinear forecasting method to extract the short-term predictability and long-term un¬ 
predictability characteristics of the gKS solution. We consider in particular two types of temporal signals, 
a local quantity, i.e. a single point in the system, and an appropriate global measure in the sense that is a 
spatially averaged quantity. We then apply several statistical and time-series methods, some of which are 
standard in chaos theory, but they have not been used to analyze the gKS solution, including the Lorenz 
return map and a nonlinear forecasting method. The latter in particular has recently been exploited in the 
context of a wide spectrum of chaotic s ystems such as expe riments in propagating flame front instabili¬ 


ties induced by buoyancy/swirl cou pling iGotoda et al 


m 


a gas-turbine model combustor [Gotoda et al. . 2011 ]. Other examples include numerical solutions of 


2ninl ] and thermoacoustic combustion instabilities 


sets of nonlinea r differential equatio ns for modeling intermittent Rayleigh-Benard convection subjected to 
magnetic force Gotoda et 120131 ] ■ and modeling flame front instabilities induced by radiative heat loss 


Gotoda et ai?.l . l2ni2bl | with emphasis on the treatment of irregular temporal time series d ata. This forecast¬ 


ing m ethod is an appropriate extension of the algorithm proposed by Sugihara and May [Sugihara fc Mavl . 
199nl ]. in the sense that it allows us to quantify the short-term predictability and long-term unpredictability 


as a function of the prediction time by taking into account the update of the library data in the phase 
space. This is crucial in order to understand the predictability properties of the system. 

We find that as the parameter 6 is increased both local and global quantities undergo a transition from 
high-dimensional chaos, characterized by scattered graphs in the phase space and a very short predictability 
time, to a periodic solution which can be predicted up to long times. At moderate values of 5 the solution 
exhibits low-dimensional chaos and is characterized by a critical predictability time that we directly relate 
to the dispersion parameter <5. We also compare the results from the gKS equation to pure stochastic 
processes like fractional Brownian motion to infer conditions under which the methodology can distinguish 
between the gKS signal from pure stochastic signals. 

Section [2] presents the different dynamical states observed in the gKS equation as the dispersion 
parameter is increased and we define there the two temporal signals mentioned earlier. Section [3] briefly 
details the nonlinear forecasting method and Sec.|3]outlines the difference between a pure stochastic process 
and a deterministic chaotic signal. The results for the local and global analysis are presented in Sec. O 
Sec. [6l and Sec. [71 We finally conclude in Sec. [8l 


2. Spatiotemporal dynamics of the gKS equation 

Gonsider the gKS equation ([2]) in a domain [—L, L], i.e. of size 2L, with periodic boundary conditions and 
discretized into Nl points. We choose random initial conditions u(x, 0) = ^{x) where ^{x) is a spatial white 
noise with {^{x)^{x')) = 26{x—x'). The numerical solution u(x, t) is obtained by adopting a pseudo-spectral 
method for the spatial derivatives that uses the Fast Fourier Transform (FFT) to transform the solution 
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Fig. 1. Spatiotemporal patterns of the gKS equation for the spatial and temporal regions of x G [—200, 200] and t G [500,1500], 
respectively, and for different values of 5. Different colours represent the magnitude of u with high and low values in yellow 
and blue, respectively. 


to Fourier space. The nonlinear terms are evaluated in real space and transformed back to Fourier space 
by using the inverse FFT. The solution is then propagated in time by making use of a modified fourth- 
order exponential time-differencing-time-steppi ng Runge-Kutta ( ETDRK4) scheme Kassam &: Trefethenl . 
2005], which is very stable for stiff equation s Shen et ai . 20111 ] . It has been reported in previous stud¬ 


ies 


Elder et al . 19971 : Wittenberg &: Holmes . 1999 ] that for sufficiently large L many modes become active 


and the transition to spatiotemporal chaos in the KS equation becomes independent of L. On this basis 
we choose L = 500 and Nl = 5000. It is also important to emphasize that by taking a large number of 
discretised points (in our case 5000) we ensure the probability density function of the random initial condi¬ 
tions converges to a Normal distribution of mean zero and variance one so that starting the computations 
with different random initial conditions has little influence on the dynamic solutions. Note that solutions 
with different means can be on th e other hand mapped to e ach other via a Galilean transformation (see 


Tseluiko Sz Kalliadasid . 20141 ] for the gKS equation and Coullet Sz Fauvd . 1995 : Malomed . 1992 ] for 


e.g. 

other systems with this type of symmetry). We also choose the value u = 1 and we vary the dispersion 
parameter 6 from 0 to 3. Note also that to improve the numerical precision, all numerical results, when 
necessary, have been obtained after averaging over 10 different realizations of the initial conditions. 

Figured] depicts spatiotemporal patterns of the gKS equation for different values of 6. Clearly, for 5 = 0, 
i.e., for the KS solution, the well-known spatiotemporal chaos is formed. As we increase 6 up to 0.2 localized 
coherent structures (solitary pulses) start to appear but with complicated chaotic dynamics due to their 
continuous interaction. Regularly-arrayed pulses with unsteady interpulse distances become predominant 
at 5 = 0.4, and for sufficiently large 6, i.e., in the dispersion-dominated regime, strongly regular structures 
with nearly co nstant interpulse distances are observed. These basic localized structures were firstly reported 


since Balmforth. 

1995; 

Chane et al. 

1995; 

Ei & Ohta. 

1994; 

Tseluiko & Kalliadasis. 

2014]. 


2.1. Definition of local and global time series 

We define here the magnitudes of interest which will be analyzed later on by means of several time-series 
tools to characterize the dynamics of the gKS solution as function of 5. We distinguish between a local 
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Fig. 2. Extraction of the temporal evolutions of u{x,t) at a; = 0 which defines the local magnitude uj^(t) for the ID gKS 
equation. Different colours represent the magnitude of u with high values (yellow) and low values (blue). 


quantity, i.e. one that represents the dynamics of a single arbitrary point of the whole domain, and a global 
quantity as a spatially averaged measure that gives information of the global dynamics. We hence consider 
the temporal evolution of the gKS solution u(x, t) located at x = 0 and define our local time signal as: 


UL{t) = u{0,t). 


( 3 ) 


Note that this is actually a standard local quantity which ha s been used f o r exa mple to characterize 
the dynamics of Turing patterns in a reaction-diffusion system Aragon et al. . 2ni2| |. Figure [2] shows the 
temporal evolution of ul for the case of 6 = 0.1. 


On the other hand, we define our global measure as the second moment of the gKS solution: 

"L 

l2. 




[u{x,t) — u{t)Ydx, 


( 4 ) 


where the overbar denotes spatial average. This quantity has been considered in other settings, 
for e xample in the nonlocal gKS e quation describing the interface evolution of conducting falling 


films [Tseluiko fc Papageorgior 


ai gKb e c 

i^ l2ninl| . 


Our study will be based on how the dispersion (controlled via 


the parameter 5) affects the dynamics of both the local ul and global uq quantities. 


3. Nonlinear forecasting method 

To study the predictability properties of the gKS equation we adopt a nonlinear forecasting method¬ 
ology ^whlch_Js_an_appropriately extended version of the algorithm proposed by Sugihara and 
May [Sugihara fc Mavl . 199fll ]. and it has been recently ap plied successfu l ly in a wide spectrum of phys¬ 
ical systems, from thermo-aco ustic combust i on inst ability Gotoda et al . 201 ll ] to flame front instability 
induced by radiative heat loss Gotoda et al . ^012bll a nd intermittent behavior of Raylegh-Benard convec¬ 
tion subjected to magnetic force [Gotoda et d. . 2013( |. The method proceeds as follows. Given a temporal 
signal, say u{t), for t G [0,T/-], we divide it into two intervals, namely t G [0,and t G {tL,Tf], corre¬ 
sponding to a library and reference set of data, respectively. The values from the library data are used to 

predict the temporal evolution of u{t) for t > tL which are then compared to the reference data^ _ 

W e construct the phase space from the time series data based on Takens’ embedding theorem F.Taken^ . 
198ll | by considering time-delayed coordinates, i.e. a vector in the phase space is given as: 


U(ti) = {u{ti),u{ti - r),..., u{ti - {D - l)r)). 


(5) 
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where i = 1,..., n with n the number of values in the data set, D is the dimension of t he phase space, and 
T is a lag time which is estimated by making use of the mutual information prop osed by Fraser &: Swinnev , 
IQSd i. in a similar way as it was also done in Refs. Gotoda et al . 2012bl. 2013l |. In this study we set D = 3 


as we find it to be the optimal value which compromises numerical convergence (i.e. results do not change 
for larger values of D) and computational time (large values of D become computationally expensive). 
We define Uj = U(tj) as the hnal point of a trajectory in the phase space, and U^, for k = 
its K neighboring vectors which are searched from all data points in the phase space. The future value 
corresponding to Ufc after some time T = NAt is denoted as u{tk + T), where At is the sampling time of 
u{t) and N an arbitrary integer. The predicted value u{tf + T) from U/ is then obtained by a nonlinearly 
weighted sum of the library data u{tk + T) given as: 


Ef=i u{tk + T)e 


u{tf + T) = 


-du 


where dk = ||U(t/) — U(tfc)||. A comparison between the predicted value u{tf + T) and the reference value 
u{tf + T) is then quantified by looking at the standard correlation factor which is defined as: 

K[u{t)u{t)] 


C = 




( 6 ) 


where K[u{t)u{t)] represents the covariance between both signals, and au and au are the standard devi¬ 
ation of u and u, resp e ctively. It sh o uld b e noted that our nonlinear forecasting method was developed 
in Refs. [Gotoda et al\ . 2011. 2012b . 2013], where the temporal evolutions of u are predicted by updat¬ 
ing the library data and so we are able to capture the determinism of the system (with keeping the size 
of the updated library data c onstant). This point is actually the difference with the Sugihara-May algo¬ 
rithm [Sngihara fc Mav . lOOfll ] where the library data is not updated. In this sense, the correlation coefficient 
C is estimated in terms of the duration tp of the actual temporal evolutions of u added to the library data 
which is something that has not been proposed before for extracting the short-term predictability and 
long-term unpredictability characteristics of chaotic dynamics. 


4. Chaotic versus stochastic behavior 


Ghaotic signals are characterized by exhibiting a high sensitivity on the initial conditions which in turn 
induces a decay of predictability in time. In particular, they exhibit short-term predictability and long¬ 
term unpredictability of the trajectories in the phase space, something which is referred to as orbital 
instability. A natural question of fundamental interest is then, given a time signal how can one distinguish 
between a pure deterministic chaotic signal from a stochastic process. Although in this work we focus on 
the deterministic gKS equation, we also analyze a stochastic signal by means of the nonlinear forecasting 
method which will be used in the subsequent sections for comparison with pure chaotic signals. 

Gonsider the fractional Brownian motion which is a Gaussian process Ba{t) with zero mean, 
i.e. E[RQ,(t)] = 0, and a time dependent variance and covariance which are given as: 


E[R„(t)R„(t)] 
E[B^{t)B^{s)] = (t“-^ 


+ s 


a—1 


- | s - t |“-')/2 


(7) 

( 8 ) 


respectively, where a is the scaling exponent (Mandelbrot fc Nessl . Il968l | . This is a well-known stochastic 
process model, the PSD of which exhibits a powe r-law behaviour of the type l//“, describing irregular 
behavior in many areas of science and engineering Biagini et al. . 2008( 1. On the other hand, the increment 
process defined as ABa = Ba{tipi) — Ba{ti) is fractional Gaussian noise. The particular case of a = 5/3 
is referred to as a Gaussian K olmogorov signal, in reference to the same type of PSD observed in fully- 


developed turbulence [see e.g. IVergassola et al . 1993 : Katnl et al. . 2001 ]. For a = 2, it corresponds to 


Brownian motion, and the increment process is white noise. 

To demonstrate the ability of the nonlinear forecasting methodology proposed in the preceding section 
for distinguishing chaos from stochastic process, we investigate the variations of the correlation coefficient 
C in terms of the predicted time tp for the two representative cases of Brownian (a = 2) and fractional 
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Fig. 3. Correlation coefficient Cftp) for (a) Brownian motion Bait) with a = 2 and the corresponding increment ABa{t) 
(b); and (c) fractional Brownian motion Bait) with a = 5/3 and the corresponding increment ABait) (d). Panels (a) and (c) 
are plots in log-log scale and panels (b) and (d) in semi-log scale. 

Brownian motion [a = 5/3) and the corresponding increment processes. We hence take u(t) = Ba{t) 
in Eq. ([6]) and compute the corresponding correlation coefficient for the two aforementioned values of a. 
The results for both cases are shown in Fig. [3j We observe that C{tp = 1) ~ 1 in both cases and it 
decreases with tp. This means that the one-step-ahead short-term prediction with high accuracy is feasible 
for Brownian/fractional Brownian motion. However, when we compute the increment process ABa the 
correlation factor is almost zero regardless of tp, indicating that it cannot be predicted. As we will show in 
the following section, this is a distinctive feature of a pure stochastic process and will allow us to distinguish 
it from the deterministic chaos of the gKS equation. 

5. Local analysis 

5.1. Transition from chaotic to periodic solutions in the gKS equation 

We start by studying how the dynamical state of the gKS solution undergoes the transition from a chaotic 
behavior, observed at low values of 5, to a pulsating periodic dynamics which appears at sufficiently large 
values of 6. To this end, we are going to perform a time-series analysis. 

5.1.1. Bifurcation diagram. 

We look at the bifurcation diagram of the gKS equation by plotting the local maxima of the temporal 
signal upf) which we define as the values upm that satisfy the following condition: 

upm = {unit) : up{t) = 0, ilpf) < 0}. (9) 

Temporal evolutions of up{t) and the corresponding bifurcation diagram for different values of 6 are shown 
in Fig. m We note that a repeated oscillation of upf) is displayed as a single point in the bifurcation 
diagram. We observe that up{t) exhibits highly irregular fluctuations at h = 0, which give rise in turn to 
the formation of uniformly distributed values of upm in the bifurcation diagram. For 5 = 0.2, up{t) starts 
to exhibit pulsating time fluctuations, and a hollow region appears in the bifurcation diagram around the 
values of 0 < upm < 2. For 5 > 0.4, these pulsating oscillations become more clear with temporal structures 
separated by irregular times, which leads to the appearance of two dominant branches in the bifurcation 
diagram. At this point the amplitude of the oscillating signal up{t) is notably amplified as 6 is increased 
which induces a gradual increase of the width separating the two branches in the bifurcation diagram. 
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Fig. 4. (a-d) Temporal evolutions of U]^{t) for different values of 5. (e) Bifurcation diagram of uj^rn- Panels (f) and (g) show 

a magnification of the selected area in panels (e) and (b), respectively. 


In the dispersion-dominated region of <5 > 2, only the upper dominant branch persists, indicating the 
formation of strongly periodic pulsed oscillations. Note however that, as shown in the top right panel of 
Fig. m the pulsed amplitude is oscillating in time which is a consequence of the weak interaction between 
pulses. For larger values of 5 such interactions would become weaker and a pure periodic signal would be 
recovered. It is important to point out that the main relevant change in the bifurcation diagram occurs at 
5 ~ 0.4 where there is an in dication of the ons e t of interior crisis . Seve ral kind of crisis phenomena in the 
KS equation were found by Chian et al . 2002 : Remoel fc: Chianl 2003| |. but there are no previous studies 


reporting the existence of such interior crisis in the temporal evolution of uiit) for the gKS equation. 


5.1.2. Power spectral density. 

We analyze next the PSD of the signal uiit) for different values of 5 which is defined as 

Suif) = {\uLim, ( 10 ) 


where uiif) represents the Fourier transform in the frequency domain of uiit) which is defined as urif) = 
{1/VT) — ul] exp {ift), where T is the time interval of observation and ul is the mean value of 

the signal UL{t) over T. The symbols (...) in Eq. ( 1101 ) denote average over different time intervals - we note 
that we are assuming statistically stationary solutions. Figure [5] shows that for 5 < 0.2 the PSD consists of 


a continuous distribution region at low frequencies followed by an exponential decay at higher frequencies. 
It is important to note th at an exponential behavio ur for the PSD is usually associated with deterministic 
chaotic systems (see e.g. Maggs fc Morales . 201 1| | and references therein) which indicates that the gKS 
solution is clearly chaotic for small values of 5. For 5 = 0.4 multiple peaks appear (2/^, 3/m, 4/m • • •) with 
basic frequency fm = 0.035 corresponding to the temporal pulsating structures of UL{t) observed in Fig. [H 
and indicate the multiple-periodicity nature of the signal. The appearance of these peaks becomes more 
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Fig. 5. Log-log plot of the PSD Su{f) of UL{t) for different values of 5. The insets in (a) and (b) correspond to a semi-log 
plot to show the exponential dependence. The dashed lines correspond to an exponential function exp(—/rg) with rg = 33 
for (5 = 0 (a) and rg = 48 for 5 — 0.2 (b). 


prominent in the dispersion-dominated region of 6 = 2. 


5.1.3. Dynamical structure of the phase space. 

We compute the Lorenz return map which is obtained by plotting pairs of successive local maximum 
points, i.e. we define the kfh. local maximum as ULm,k and we plot ULrn,k+i against ULrn,k- Figure Ha) 
depicts the results for <5 = 0, 0.2, 0.4 and 2. An entirely scattered structure appears for <5 = 0, indicating 
the presence of chaotic behavior. For 6 = 0.2 we get a map which is mainly gathered in three or four 
regions with scattered structures, which with increasing 6 up to 0.4, they progressively group in four main 
regions with locally scattered plots. Interestingly, two of them follow the straight line ULm,k+i = ULrn,k- 
In the dispersion-dominated regime with 5 = 2, these four points aggregate at one region which belongs 
to the line ULrn,k+i = ULm,k- This is evidence of the formation of strongly periodic pulsating oscillations, 
but note that the group located at the top right of the map is still exhibiting a scattered structure. Note 
also that the amplitude of ULm in Fig. 0 ] yields the entirely scattering plots in the localized region. This 
indicates that the dynamical state of tt at 5 = 2 is dominated by a strongly periodic pulsed oscillation, but 
it involves small fluctuations in the oscillation amplitude (see top right panel of Fig. 0|). 

We also look at the three dimensional (3D) phase space which we construct by considering 
[uL{t),UL{t — T),UL{t — 2t)]. The results for different values of 5 are plotted in Fig. Hb)- As before, we 
observe a entirely scattered graph for 5 = 0 which fills the core of the attractor, indicating hence the 
presence of high-dimensional chaos. It still remains a scattered structure for 5 = 0.2, but is different from 
that at 5 = 0. The difference between the two cases is hard to be extracted by the conventional linear 
analysis such as the PSD (see Fig. 4). In this sense, drawings of the Lorenz return map and the 3D phase 
space are useful for capturing the subtle differences between the dynamical state of the KS and the gKS 
equation with small 5. With 5 increasing up to 0.4 the overall shape of the attractor exhibits the geomet¬ 
rical structure of a strange attractor which is a consequence of the occurrence of multiple-periodic pulsed 
oscillations with irregular time intervals, indicating hence the presence of low-dimensional chaos. For 5 = 2, 
we observe the emergence of a limit cycle with some small fluctuations corresponding to the fluctuating 
amplitude (see Fig. [4]). 
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Fig. 6. (a) Lorenz return map of U£(t) for different values of (b) 3D phase space lu£(t),ui^(t — T),ui^{t — 2r)] for <5 = 0 

(r = 2),S = 0.2 (r = 2),S = 0.4 (r = 6), and <5 = 2 (r = 7). 


5.2. Nonlinear forecasting analysis 

We apply the methodology described in Sec. [3] to study the predictability properties of UL{t). We divide 
the time series UL{t) into library and reference data by considering the values of ul up to = 2500 to be 
used as library data and those thereafter as reference data. 

Figure El^a) depicts the correlation coefficient C as function of the prediction time tp for different 
values of 5. We observe that for the KS equation ((5 = 0), C at tp = 1 is nearly unity, which means 
that one-step-ahead prediction of the dynamical state is feasible with high accuracy. However, C decreases 
rapidly with tp reaching nearly zero at tp = 100 . As 5 is increased, we observe that a short plateau of 
constant C starts to appear which, after a critical predicted time, say tp^c, is followed by a rapid decay the 
rate of which depends on the particular value of 5, and which can be approximated by a power-law tail, 
C{tp) ~ t~p, with the exponent 7 depending on 6. In particular we hnd that 


C(tp) ~ 


1 if tp < tp^c) 
if tp > tp^c, 


( 11 ) 


with 7 = 0.67±0.06 for 5 = 0, 7 = 0.55±0.07 for <5 = 0.1, and 7 = 0.41±0.04 for 5 > 0.2. Based on the above 
relationship, the predicted critical time tpp is hence numerically obtained by htting a power law at long 
times and a constant function at short times, and find the point where these two behaviours meet [see the 
inlet panel of Fig. [71[b) for an example]. The fact that the same exponent 7 is observed for 6 > 0.2 suggests 
data of C for different values of 6 should collapse when plotted against the ratio tp/tp^c- The results are 
presented in Fig.[7][b), where we can observe an excellent data collapse indicating indeed that for 6 > 0.2 the 
exponent 7 becomes independent of 6. Note that the appearance of such constant region followed by a rapid 
decrease is the manifestation of short-term predictability and long-term unpredictability characteristics of 
chaos which are separated by the critical value tpp. It is important to emphasize that even in the case of 
6 = 0.4, where the spatial solution is characterized by localized structures, the temporal behavior exhibits 
a chaotic motion which is a consequence of the irregular interaction between these structures. This regime 
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Fig. 7. (a) Log-log plot of the time correlation coefficient C{tp) for up{t) against the predicted time tp for different values 

of 5. The inlet panel corresponds to the correlation coefficient for a pure periodic function sin(27r/t). (b) Data collapse of 
C against tp/tp (. where tp ,, is the critical predicted time after which C starts to decay. The dashed line corresponds to a 
decaying power-law tail with exponent —0.4 ± 0.03. The inlet panel shows how the critical predicted time tp^^ is defined. 


corresponds to low-dime nsional chaos, or what sometimes is referred to as weak/dissipative turbulence (in 
the “Manneville sense” Mannevillel . 1985], as noted in the Introduction). 

In the dispersion-dominated region with <5 = 2, the region of constant C becomes significantly larger 
followed by a slower decay. We note that if the dynamical state is periodic oscillations, C is unity regardless 
of tp meaning a long-term predictability (see inlet panel of Fig. [7Ka) for an example of temporal prediction 
of a sine wave sin(27r/t) with / = 0.035). The gradual decay of the gradient of C at d = 2 is attributed 
to small fluctuations in pulsed oscillation amplitude extracted by the Lorenz return map, i.e. uniformly- 
scattered plots in a single region. 

We also scrutinize the increment process of up which we define as Aup = up(tj+i) — up(tj). Figure 
[8] shows the correlation coefficient C for 5 = 0, 0.2 and 0.4, where we observe a similar trend as we did 
with the results for up [cf. Fig. [TKa)]. We note that the behavior of Aup is completely different to the one 
observed for the increment process of the fractional Brownian dynamics (see Fig. [3|). This indicates that 
a way of distinguishing between chaos and pure stochasticity is by looking at the predictability properties 
of the increment process, as opposed of looking at the original signal - this is particularly true for high¬ 
dimensional chaos (5 = 0) as both signals exhibit a similar trend for C{tp). 

Based on these results, the short-term predictability and long-term unpredictability of the gKS equation 
are summarized in Fig. [9l As we increase 5, the dynamical state is stabilized and the gradient of the 
correlation coefficient against the prediction time, i.e. the exponent 7 , decreases. The degree of short-term 
predictability can then be roughly estimated by the critical prediction time tp^c- Figure fTOl shows this time 
and the exponent 7 as function of 6, observing that tp^c gradually increases for 6 > 0.2 and 7 rapidly 
decreases for 6 < 0.2 and remains constant around the value 0.4 until d ~ 1, after which starts slowly 
decreasing. It should be pointed out however that for these large values of 6, the power-law behaviour at 
long times (and in turn the exponent 7 ) is largely affected by cross-over effects: to find the truly power-law 
decay we would need to go to longer values of tp. These three regimes are identified with the existence 
of high-dimensional spatiotemporal chaos at low values of 6, low-dimensional chaos at moderate values, 
and periodic solutions at large values. These results demonstrate that the proposed nonlinear forecasting 
methodology allows us to fully characterize the dynamics of the system in terms of its predictability 
properties. 
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Fig. 8. Correlation coefficient C{tp) for the increment process Aup{t) and for different S. 


Short-term predictability Long-term unpredictability 



Fig. 9. Schematic representation of short-term predictability and long-term unpredictability characteristics. 


6 . Global analysis 

We repeat the same analysis as in Sec. [5] but considering the global magnitude ucit) as defined via the 
second moment in Eq. Q. Figure [TTlal depicts the temporal evolutions of ucit) for different 6. Irregular 
fluctuations are clearly formed for 6 < 0.2 and we now begin to observe some pulsating oscillations with 
irregular time interval separations at d = 0.4 which become more regular at <5 = 2. It is important to 
emphasize that in contrast to what we observed for the local magnitude up [cf. Fig. 0], these pulsating 
oscillations are superimposed with a long wave oscillation of irregular frequency. 

We compute the PSD of ucit) observing a similar behavior as for the local magnitude up when 6 is 
increased [see Fig. fTTIb)]. We also look at the Lorenz return map and the 3D phase space observing that 
now the transition from a completely scattered graph to structured graphs appears to be at larger 6 values 
(see Fig. d. Figure fT3l shows the correlation coefficient C for the signal uq and for the increment process 
Aug = uciti+i) — uciti), respectively. As before, we observe an initial regime of short-term predictability 
followed by a power-law decay of C indicating the long-term unpredictability which is decreased as we 
increase 5. 

We note that the formation of periodic pulsed oscillations which are now superimposed on a long wave 
of irregular frequency (see Fig. [TT] for 5 = 0.4) makes the signal uq more chaotic as compared to up, in 
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Fig. 10. Variations in and 7 = |dlnC/dlntp| for uj^it) as a function of 5. Typical errors are 3 and 0.03 for tp,, and 7 , 
respectively. 



Fig. 11. (a) Temporal evolutions of UQ{t) for different <5. (b) PSD of the global magnitude uoit) for different 5. 


the sense that now the transition to long-term unpredictability appears before than in the case of ml. This 
is more clear if we plot the critical predicted time tp^c against 6, observing that now it increases until 
5 ~ 1.1 after which it remains constant until 6 = 2.0 (cf. Fig. fT^ . As for the gradient |dlnC/dlntp|, we 
observe that it initially decreases for 6 > 0.2 after which it remains constant until 6^2 for which it starts 
decreasing again. It is noteworthy that the critical predicted time is shorter than that for the local signal 
indicating hence that the global quantity appears to be less predictable. 


7. Permutation entropy 

Finally, we present an alternative analysis of both the local and global signals, up and uq, respectively, 
based on the concept of permutation entropy, which is an invariant measure of the complexity of the 
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dynamics Bandt fc Pompe . 2002j |. and has been shown to be a good candidate to e.g. quantifying the 
dynamics in combustion instability [Gotoda et ai . 2012a |. This analysis complements the previous results 
and provides more information about the difference between the magnitudes of the local and global sig¬ 
nals. A m odified vers i on of the permutation entropy has also been recently proposed to study heartbeat 


dynamics Bian et ali . I2ni2l | 


The permutation entropy hp can be estimated by quantifying the degree of randomness from a sequence 
of ranks in the values of a time series. In particular, given a sequence with embedding dimension Df., we 
index all possible permutations (T*eO of order > 2 as tt. Each permutation represents a coarse-grained 
pattern of the dynamics when the sequence consists of successive data points taken from a time series. We 
hrst count the number of permutations denoted as ^(vr) for all vectors X(t) = {u{t), u{t+l ),... , u{t+De — l)) 
consisting of sequences of order D^. We then calculate the relative frequency for each permutation to obtain 
PeiT^) = q{Tr)/{N — {De + 1)). We Set De = 7 in this study, and following the definition of Shannon entropy, 
hp is then calculated as: 


hp — 


-EPe(^)log2Pe(vr) 

TT 

log2 Del 


( 12 ) 


where log 2 T*e! corresponds to the maximum permutation entropy and is used for normalization so that 
0 < /ip < 1. Note that with this definition the lower bound {hp = 0) corresponds to a monotonically 
increasing or decreasing process, while the upper bound {hp = 1) corresponds to a completely random 
process. 

Figure 19 shows the results of hp for both and uc as we change the values of 5 from the high¬ 
dimensional chaos region {6 = 0) to the low-dimensional chaos region (5 = 0.4). We observe that hp 
decreases for both ul and uq with increasing 5, as a consequence of the dispersion mechanism that arrests 
chaotic behavior. It is important to note that hp is significantly larger for the case of ug, indicating that 
the dynamics of the global measure is more difficult to analyse than the corresponding local one, something 
which is a consequence of the averaging procedure. 


8 . Conclusions 

We have presented a detailed and systematic study of the dynamics of the spatiotemporal behavior of 
the gKS equation by making use of statistical methods and time-series analysis based on chaos theory, 
including a nonlinear forecasting method. The gKS equation is a generalization of the KS equation to 
include dispersion, and has been used as a prototype for the study of pattern formation dynamics and 
spatiotemporal complexity in active dispersive-dissipative media. The spatiotemporal solution of the gKS 
equation exhibits a rich dynamics that strongly depends on the parameter 5 controlling the strength 
of the dispersion. In particular, one observes spatiotemporal chaos at low values of 5, the emergence of 
localized structures (pulses) constantly interacting with each other at moderate values (regime also known 
as interfacial turbulence), and periodic solutions at large values of 5. 

We have defined two different temporal magnitudes to be analyzed; a local signal, the extracted solution 
in the middle of the system, and a global measure, the second moment of the solution. This allows for 
two different ways, a local and a global, of accessing the dynamical properties of the system providing 
then a complete picture (usually in experiments one has access to one of the two only). The local analysis 
has revealed that the high-dimensional chaos regime (5 < 0.2) is characterized by a completely scattered 
graph in the phase space, a power-law decay PSD, and a very short predictability time. As the dispersion 
parameter is increased, the spatiotemporal solution starts to exhibit the emergence of solitary waves (pulses) 
which continuously and chaotically interact with each other. This corresponds to a low-dimensional chaos 
regime exhibiting a structured graph with a strange attractor in the phase space, and with a typical 
critical predictability time up to which the solution can be predicted. For larger values of 6 the critical 
predictability time starts to increase as the solution approaches a periodic steady state {6 > 1.8). 

We observe similar behavior for the global analysis but the transitions from high-dimensional to low¬ 
dimensional chaos and finally to periodic solutions occur at larger values of 6, namely at 5 ~ 1 and d ~ 2.0, 
respectively. This is a consequence of the global spatial average we consider since any sign of localized 
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Fig. 12. (a) Lorenz return map of UQ{t) for different 5. (b) 3D phase space of the global magnitude UQ{t) for <5 = 0 (r = 2), 

S = 0.2 (r = 2), 5 = 0.4 (r = 3), and 5 = 2 (t = 8). 




Fig. 13. (a) Correlation coefficient C(tp) of uo{t) (a) and AuQ{t) (b) for different values of 5. 

structure is lost during such averaging procedure. In this sense we can conclude that the global measure 
is more chaotic than the local signal. This is also confirmed by looking at the permutation entropy of the 
temporal signals observing that the global signal has larger entropy than the local one. 

We have also contrasted the results from the deterministic gKS equation to a pure stochastic process 
(fractional Brownian), observing that one needs to look at the increment process of the signal in order to 
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5 

Fig. 14. Variations of tp ,. and 7 = |dln C/d lntp| for UQ{t) as a function of 5. Typical errors are 4 and 0.03 for tp^c a^nd 7 , 
respectively. 



Fig. 15. Variations of hp for up and uc computed by using different values of J. 


distinguish between stochastic and high-dimensional chaotic signals. In particular, we have observed that 
while both the gKS solution at low values of 6 and the fractional Brownian motion exhibit a very short 
predictability time for the original signal, the increment process shows zero predictability properties for the 
stochastic process and short predictability properties for the chaotic signal. To the best of our knowledge 
this is the first study to provide a comprehensive interpretation of the dynamical state of the gKS equation 
from the viewpoint of nonlinear forecasting. 

Finally, there is a number of interesting questions related to the analysis presented here. For instance, 
how can one use elements from stochastic processes (e.g. in Krumscheid et all 2ni.‘ll ]) to appropriately 
modify the nonlinear forecasting method when the gKS equation is postulated for a nonlinear process for 
which the precise underlying model is not available. Another interesting problem would be the extension 
of the methodology proposed here to more involved equations, such as those describing the dynamics of 
the falling film away from criticality. We shall examine these and related issues in future studies. 
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